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Abstract. Shock waves in molecular clouds should evolve 
into continuous or C-type structures due to the magnetic 
field and ion-neutral friction. We here determine whether 
and how this is achieved through plane-parallel numeri- 
cal simulations using an extended version of ZEUS. We 
first describe and test the adapted code against analytical 
results, laying the necessary foundations for subsequent 
works on supersonic ambipolar diffusion, including C-type 
jets and shock instability. 

The evolution away from jump shocks toward the nu- 
merous steady C-shock sub-types is then investigated. 
The evolution passes through four stages, which possess 
distinctive observational properties. The time scales and 
length scales cover broad ranges. Specific results arc in- 
cluded for shock types including switch, absorber, neu- 
tralised, oblique, transverse and intermediate. Only in- 
termediate Type II shocks and 'slow shocks', including 
switch-off shocks, remain as J-type under the low ion levels 
assumed. Other shocks transform via a steadily growing 
neutral precursor to a diminishing jump. For neutralised 
shocks, this takes the form of an extended long-lived ramp. 

Molecular hydrogen emission signatures are presented. 
After the jump speed has dropped to under 25kms _1 , a 
non-dissociative jump section can dominate the spectra for 
a long period. This produces a high-excitation spectrum. 
Once the jump has further weakened, to <8kms _1 , the 
fully developed ion front is responsible for brisk progress 
towards a constant C-type excitation. The time scale for 
emission-line variations is ~ (6/rii) yr, where n^ is the pre- 
shock ion number density. 
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1. Introduction 



C-type shocks are frequently invoked to interpret the sig- 
natures of accelerated and excited gas in molecular clouds. 
These magnetically-mediated, two-fluid shocks are able to 
explain the survival of the molecules in shocks with ve- 
locities exceeding 20kms~ 1 as well as the low excitation 
of observed spectra (e.g. the strong atomic fine-structure 
lines of oxygen and the low values of the 2-1 S(l)/l-0 
S(l) ratio of molecular hydrogen). However, several basic 
questions on the existence of C-shocks remain to be an- 
swered. Under what conditions do steady-state C-shocks 
form? Is there sufficient time for the various C-shock con- 
figurations to be realized? Alternatively, how can we rec- 
ognize an incomplete or proto-C-shock? Here, as part of a 
wider program to investigate supersonic ambipolar diffu- 
sion, we develop a numerical model and apply it to study 
the time-dependence of planar C-shock configurations. 

A theory for steady C-shocks in molecular gas was con- 
structed by Draine (1980) and applied by Draine, Roberge 
& Dalgarno (1983). The combination of low ion fraction, 
strong cooling, and significant (but not high) magnetic 
field results in a shock in which ion-neutral drag provides 
the viscosity, cooling keeps the gas supersonic and the field 
via the ions provides an extended cushioning layer which 
inhibits molecular dissociation (Draine & McKee 1993). 
Steady solutions for transverse- field C-shocks were further 
analysed by Chcrnoff (1987), Roberge & Draine (1990) 
and Smith & Brand (1990a). Smith, Brand & Moorhouse 
(1991) looked at the high-field 'shock absorbers'. Wardle 
& Draine (1987) & Smith (1993a, b) presented a theory 
for oblique-field C-shocks, including the parallel-field C- 
switches. Flower, Pineau des Forets & Hartquist (1985) 
considered the chemical aspects of steady C-shocks, work 
which has developed into a theory for C-shock chemistry 
(see Flower et al 1996). Predictions for steady, planar C- 
shocks were presented, in addition, by Smith & Brand 
(1990b), Smith (1991, 1995) and Kaufman & Neufeld 
(1996a,b). 
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Are these steady-state predictions relevant? There are 
two other possibilities: instability and evolution. Insta- 
bility has been investigated by Wardlc (1990, 1991a,b), 
Toth (1994,1995) and now by Stone (1997) and Mac Low 
& Smith (1997a, b). Conditions for stability are either (i) 
neutral Alfven numbers M a — v s /v a < 5 ,where v s and v a 
= B /(47T / 9„) 1 / 2 are the shock and Alfven speeds, respec- 
tively, or (ii) rapid recombination so that the ion fraction 
is fixed locally rather than via advection. However, it is 
not clear that the steady state is approached even in the 
absence of instability. This evolutionary question, which 
has not been examined before (Pineau des Forets 1997), is 
the focus of this paper. Furthermore, we predict the speed 
and character of changes in emission line strengths that 
identify a proto-C-shock, and distinguish these changes 
from those produced in an unstable C-shock (as evalu- 
ated by Neufeld & Stone (1997) and Mac Low & Smith 
(1997b)). 

Our attempt to evolve flow patterns from a sharply dis- 
continuous J-shock to a C-shock under constant ionization 
conditions may not exactly correspond to a physical situa- 
tion. Rather, we use these examples to envisage how a flow 
accomodates to changing conditions such as those caused 
by the start of a stellar outflow, the impact of a jet on a 
cloud, or a cloud-cloud collision. Our primary goal here is 
to obtain a deeper understanding of how molecular shocks 
behave dynamically. We do use an alternative initial con- 
dition, with a smooth transition following a hyperbolic 
tangent function, to test the dependence of our results on 
the assumption of an initially sharp discontinuity. 

This paper also represents one step in our exploration 
of supersonic ambipolar diffusion. The original ZEUS am- 
bipolar diffusion code of Mac Low et al (1995) has been 
extended to cover more general molecular cloud condi- 
tions, by including ion mass conservation as opposed to the 
ealier assumption of a fixed ion number density. Hence it 
is first necessary to test the behavior of the extended code 
against standard C-shock solutions (§2). The fluid remains 
as before isothermal. We thus concentrate on the dynam- 
ical aspects. The 'cool C-shock approximation' (Smith & 
Brand 1990a) allows us to extract quantitative predictions 
for temperature and line emission. 

2. Framework 

2.1. Numerical Methods 

For our numerical computations, we use a modified ver- 
sion of the ZEUS code 1 (Stone & Norman 1992a,b). Am- 
bipolar diffusion was added to ZEUS by Mac Low et al. 
(1995), who described the basic interface with the ZEUS 
code. Summarizing, that work made four approximations: 
isothermality of ions, electrons and neutrals, no electron- 
ion drift, ion density dependent in power-law fashion on 

1 Available for community use by registration with the Lab- 
oratory for Computational Astrophysics at lca@ncsa.uiuc.edu 



neutral density, and no ion inertia or pressure. This al- 
lowed us to neglect, respectively, the energy equation, 
Ohmic diffusion, and the equations of ion mass and mo- 
mentum conservation. This approach has proved adequate 
for modelling protostellar disks in the absence of strong 
shocks. However, the flow of ions can be important in 
C-shocks, for example in the Wardle instability, which is 
driven by the flow of ions along buckling field lines in the 
shock front. 
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Fig. 1. The standard transverse-field C-shock reaches this 
steady state after 5 10 10 s. The grid resolution R= Ax/L n , 
where L„ is the neutral deceleration length, roughly corre- 
sponding to the shock thickness. Here we show a low-resolution 
run with R=13 to emphasize the accuracy of the code. The dot- 
ted lines are the corresponding analytical solutions. The total 
compression S = 34.86. 



Neglect of ion inertia and pressure had allowed the 
replacement of the ion momentum conservation equation 
by an algebraic equation expressing the balance between 
Lorentz forces and ion-neutral drag in determining the 
drift velocity between ions and neutrals. This approach 
is physically accurate and allows time steps determined 
by the equivalent of the Courant condition for ambipolar 
diffusion, At < ir 1Pl p n {Ax) 2 /|B| 2 (Mac Low et al. 1995). 
Both Toth (1994) and we have found that this approach 
can be numerically unstable in the presence of steep ve- 
locity gradients as occur in C-shocks. However, in the one- 
dimensional models shown in this paper, this only occurs 
for switch shocks, and even there is not deadly to the com- 
putation, so we defer addressing this issue to our report 
on multi-dimensional models (Mac Low & Smith 1997b). 

We treat the ions as a separate fluid in the code, us- 
ing the standard ZEUS algorithms to update them. The 
equations we solve in the current version of the code are 
then the neutral and ion continuity equations, the neutral 
momentum equation, and the induction equation, as well 



Smith & Mac Low 



3 



as satisfying the zero-divergence criterion: 
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where the subscripts i and n refer to the ions and neu- 
trals, p, v, and P are density, velocity and pressure for 
each fluid, B is the magnetic field, and 7 is the collisional 
coupling constant between the ions and neutrals. 

We used the analytic jump conditions for a magne- 
tized shock (e.g. Priest 1982) to set the values of the flow 
variables on both boundaries needed to hold the shock sta- 
tionary on our grid. These are traditionally called 'inflow' 
boundary conditions, although on the downstream edge of 
our grid the flow is out rather than in. 

2.2. Analytical background and parameter definitions 

We begin by setting up a one-dimensional J-shock. We 
generated this solution from the standard shock jump 
conditions with magnetic field (e.g. Priest 1982). This is 
equivalent to assuming that the ion fraction was initially 
high, forcing the ions and neutrals to move at identical ve- 
locity. Given the shock speed, Alfven speed, temperature 
and field direction, we can define the downstream values 
in terms of the upstream values. The supersonic inflow 
and subsonic outflow boundary conditions are then set so 
that the shock front remains stationary on the computa- 
tional grid. Note that the shock front thus begins as a 
truly discontinuous jump rather than being spread out by 
numerical viscosity. We also use an initial condition with 
a smooth transition proportional to tanh (x/L) for a given 
L to check the importance of the initial discontinuity. 



Table 1. Initial and boundary conditions: the standard set. 



We start with a standard set of initial conditions. 
This standard set represents a high Alfven number shock 
(M a = 25) with a transverse magnetic field. In Table I we 
present all the physical conditions of the models shown in 
Figures 1-4, apart from the sound speed of 100 cm s _1 , 
helium abundance of 0.1, ion mass m t = 10m p , neutral 
mass m„ = 7m p /3 (assuming a fully molecular gas), and 
the coupling constant 7 = < aw > /{rrii + m n ) — 9.21 
10 13 cm 3 s _1 gm _1 , corresponding to a fixed momentum 
transfer rate coefficient of < aw >— 1.9 10~ 9 cm 3 s" 1 (a 
rather uncertain parameter). 

We then determine the corresponding steady-state 
C-shock solution. We describe this with the variables 
r(x,t) — v n /v s , q(x,t) = Vi/v s and the fixed tempera- 
ture r = kT n /(m n v 2 ). Hence r = q = 1 upstream and the 
neutral and ion compression ratios are 1/r and 1/q, re- 
spectively. Assuming ion conservation and isothermality, 
the steady-state C-shock solution is then fully described 
by the total momentum equation 
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(Smith & Brand 1990a) where L n = v a /(jpi) is the neu- 
tral deceleration scale. For the 'neutralised' isothermal C- 
shock, in which the ion abundance is a constant, the mo- 
mentum conservation condition is unaltered and the drag 
relation is then 
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a This is the total hydrogen nuclei number density. 
b The initial upstream Alfven speed is 2 kms" 1 . 



c The number of ions relative to uh- 



Full steady-state solutions are derived by numerically solv- 
ing these simple differential equations to any desired ac- 
curacy. 

2.3. Scaling 

Besides the length scale L n , we define the ion deceleration 
scale Li = L n /M a . A further scale defines the ion- neutral 
interaction region. In cool shocks, the maximum stream- 
ing speed occurs for q = (on putting t = in 
Equation 6). The value of the maximum is (r — q) m ax = 
l+M- 2 -3Ma 2/3 /2. It follows that high streaming speeds 

1 /3 

are maintained over a length scale L str ~ L n /Ma'°. Fi- 
nally note that the length scale defined by Wardle (1990, 
1991a, b) is L 8 hk = V%L n . We define here what proves 
to be a better measure of the shock length for transverse 
shocks: L sm = L n /(r — q) max . This enables low M a shocks 
to be also included 

The time-dependent shock flow pattern is fully deter- 
mined by four control parameters: the Alfven number M a , 
the Mach number A'/, the ion fraction and the initial field 
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orientation. These are the only parameters which remain 
after scaling equations (1-5) to the length scale L n and the 
time scale tfi ow = L n /v a . Therefore only one parameter 
is necessary to describe the whole class of transverse, low- 
ionisation, cold flows: a single flow simulation with a fixed 
Alfven number is relevant to a wide range of conditions. 

One must, however, ensure that the thermal pressure 
gradients remain small. This is easily achievable when 
close to the steady state since then the pressure gradient 
term r/r = c 2 s /(v n v s ) is limited by the maximum com- 
pression (S ~ V2M a for high M a ) to the value (c 2 Jvl)/S, 
where c s is the sound speed. This is expected to be well 
below unity (and therefore ignorable) in strong shocks 
in molecular clouds. However, in time-dependent flows 
in general, depending on the imposed conditions, neutral 
speeds may approach zero and even an isothermal flow 
may possess high pressure gradients. The present simu- 
lations then require holding both the Mach and Alfven 
numbers fixed. 

Furthermore, the freedom to scale parameters while 
holding M a fixed becomes physically invalid if the ion 
fraction depends on the other shock parameters. This lim- 
its the scaling regime to flows with ion-neutral streaming 
speeds less than ~ 42 kms -1 to avoid runaway ionisation 
(Draine et al 1983, Smith & Brand 1990a). 

The flow patterns presented below can thus be consid- 
ered quite general provided the above rules are not vio- 
lated. No scaling was attempted, however, when calculat- 
ing the emission line properties, in order to avoid possible 
mis-matching when piecing together the results from sep- 
arate models for the emission from the J-shock and the 
continuous section. 
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2.4- Test of accuracy 

In our time-dependent computations, J-shocks indeed 
evolve to the analytical, steady C-shock solution to within 
our numerical accuracy. Figure 1 presents the flow pa- 
rameters at a time 510 10 s. The length scale is L n = 
1.310 15 cm and the flow time scale is tfi ow = L n /v a = 
6.5 10 9 s. We find that it takes several flow time scales be- 
fore the flow pattern finally settles down to the analytical 
steady state solution. 

Figure 2 demonstrates the dependence of the time- 
dependent solutions on the grid resolution. Even a grid 
size of Kg = 50, with the number of zones across the shock 
R = Ax/L n = 13, is sufficient for dynamical purposes with 
the major error occurring in the final deceleration zone (at 
about 3 10 15 cm). The error functions analysed in Figure 2 
are the momentum error (solid lines) given by 



£m(*£) 



r+l/(2M a V 



1 + 1/(2MJ) 
and the streaming error (dotted lines) given by 



e s {x) 



r + q 



(9) 



(10) 



Fig. 2. The errors in the final configurations for simula- 
tions of the standard C-shock at the four indicated resolutions 
R — Ax/L n . The number of zones in the full grid is given by 
R g . The solid line is the momentum error and the dotted line 
is the streaming error. 



where [r s — q s ] is the analytical steady streaming solution. 
The expected second order convergence (Mac Low et al 
1995) occurs. 

3. The Evolution from J-shock to C-shock 

3.1. The four stages for transverse shocks 

We identify four stages in the evolution from a J-shock to 
a C-shock. Each stage lasts approximately 10 times longer 
than the previous one. 

— Stage 1 involves rapid ion motions. It is characterised 
by high-speed, ion-magnetosonic wave motions over a 
short time. During the first 0.002i/; OMJ <~ 10 7 s, the 
precursor moves upstream at a speed of order 1000 
kms -1 , the ion-magnetosonic speed. As shown in Fig- 
ure 3, the ions move upstream (negative velocities), 
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stoge 1: rapid ion motions 
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stage 2: ion-precursor formation 




distance (10 



Fig. 3. The initial phase in the J to C-shock transition involves 
rapid ion motions. The neutral (solid line) and ion (dash-dotted 
line) velocities are shown. Note the distance scale has been 
expanded, so that the neutral jump shock, spread over 2-3 
zones by numerical viscosity, appears rather broad. 



Fig. 4. The second phase in the J to C-shock transition. The 
ion and neutral fronts separate at speeds in the range 10-100 
km -1 . The neutral (solid line) and ion (dash-dotted line) ve- 
locities are shown, as well as 0.001e m (dotted line), where the 
momentum error e m , given by equation 9, tracks the deviation 
from the steady-state solution. 



away from the original shock front, with a speed well 
above the shock speed. The ion compression wave is 
separated from the neutral front by an ion expansion 
wave, as expected in this shock-tube experiment. 

— During Stage 2 the ion front and the neutral front 
separate over a time of 0.03tfi ow ~ 210 8 s, at inter- 
mediate speeds (Figure 4) . The ion precursor develops 
fully. The neutrals possess a weak precursor but are 
still hardly altered. 

— Stage 3 represents the approach to a C-type structure, 
as the neutrals evolve towards a fully continuous flow. 
The changes now occur at the Alfven speed of 2 km s _1 
over a time of 0.6tfi ow <~ 4 10 9 s. 

— During Stage 4 only small changes to the flow pattern 
occur. It takes several Alfven speed crossing times be- 
fore the steady state solution is finally reached, as de- 



scribed in the previous section (see Fig. 1). The final 
length scale is <~ 1.610 15 cm, approximately equal to 

3.2. Low Alfven number shocks 

Strong magnetic fields cushion a shock by broadening the 
transition region and reducing the drag heating in the 
shock front. This is especially important for C-shocks since 
the cushioning changes the character of the whole transi- 
tion by reducing the ion-neutral streaming speed whereas 
for J-shocks it is usually only felt in the downstream 
compressed gas. Hence we call these C-shocks 'shock ab- 
sorbers' (Smith et al 1991). This cushioning has observa- 
tional implications: molecules arc not so easily destroyed 
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stoge 3: neutral tronsformotion high-field obsorber 
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Fig. 5. The neutral jump shock finally disappears and the 
shock width further increases in Stage 3. The neutral (solid 
line) and ion (dash-dotted line) velocities are shown, as well as 
0.1e m (dotted line), where the momentum error e m is given by 
equation 9. 



Fig. 6. At low Alfven numbers, the ions and neutrals gently 
stream against each other. In this model, the shock speed v s — 
5 kms -1 , the neutral Alfven number M a =2.5, the Mach num- 
ber M = 50, and the compression S = 3.07. Note the distance 
scale has been doubled. The grid size R g = 200, so the reso- 
lution R = 110 zones in L n . The neutral (solid line) and ion 
(dash-dotted line) velocities are shown, as well as 0.1e m (dotted 
line), where the momentum error e m is given by equation 9. 



by ion collisions because the ion-neutral streaming speed 
is a relatively small fraction of the shock speed. In order 
to model a shock absorber, we simply reduce the shock 
speed from the standard case by a factor of 10 to simulate 
a M a = 2.5 shock. 

We find that all the transition stages are identifiable 
but last several times longer (Figure 6). Indeed it takes 
<~ 10 4 years to reach the steady state, comparable to the 
age of many outflows. Note that the total compression is 
only about 3 in the example shown. The final total length 
of the shock is ~ 4 10 15 cm. It is clear that L sm is indeed 
an accurate measure of the shock length scale. 



3.3. Oblique field 

Oblique shocks develop from the initial state to a steady- 
state, oblique C-shock without any surprises. By oblique 
shock, we specifically mean that the magnetic field is in- 
clined to the shock velocity, but not so far inclined that in- 
termediate shock solutions become possible (Smith 1993c; 
see below). In this regime, the ion and neutral speeds 
along the flow direction, qv s and rv s , behave as before, 
as shown in Fig. 7. The ion transverse speed qj,v s exceeds 
the neutral transverse speed r y v s . Note also that the ion 
transverse speed passes through a maximum, and that the 
shock is narrower than the transverse solution. Thus the 
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grid size was reduced to preserve approximately the same 
numerical resolution. 



is possible that in complex, multi-dimensional, shock con- 
figurations a mixture of the three possible shock solutions 
can be maintained at field angles less than ~ l/M a ra- 
dians, the switch- type solution is the only evolutionary 
solution in plane-parallel flows (Kennel 1988). 



oblique shock tronsformotion 
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Fig. 7. An oblique shock with the upstream magnetic field 
at an angle of 30° to the shock velocity. The ion (dash-dotted 
line) and neutral (solid line) velocity components along the 
flow direction, qv s and rv s , are shown as before, as well as the 
transverse speeds q^Vs (dashed line) and r^v,, (dotted line). 
The upstream transverse speeds were set to zero. Here R 9 = 
200, and R = 109. The shock speed is 20kms~ 1 , Mach number 
M — 200 and M a — 10, giving a total compression S = 23.1. 



switch shock 
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Fig. 8. The evolution of a switch shock structure with shock 
velocity v s — 10 kms~ , magnetic field at an angle of 2° to v s , 
a compression of S = 22.31, and M a = 5. The parallel field is 
constant (3.43 10~ 4 G) while the transverse field is 'switched 
on' from 1.20 10~ 6 G to 2.37 10~ 3 G. The ion and neutral ve- 
locity components are shown as described in Fig. 7. The grid 
size is R s = 300 zones, giving a resolution R = 121. 



3.4- Near-parallel field 

When the magnetic field is quasi-parallel, the flow does 
not reduce to the hydrodynamical equivalent, as is often 
assumed. Even as the transverse field approaches zero, the 
flow should approach a switch shock solution. This is due 
to the strong cooling in molecular shocks, or, here, the 
assumption of isothermality (Smith 1993a, c). Although it 



To compute switch shocks we had to solve several prob- 
lems. First, we had to set up an appropriate grid. If we 
take a magnetic field angle of 2° to the flow direction, 
our standard parameters lead to resolution problems: the 
neutral shock width is reduced to <~ L n /M a , but, in con- 
trast, the ion precursor extends far upstream. We solve 
this problem simply by extending the grid upstream using 
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a ratioed grid, and increasing the total number of zones. 
This problem can still be seen in Fig. 8, where we have re- 
duced M a to 5 (v s = lOkms -1 ). The whole ion precursor, 
as evident in the upstream transverse ion speed, is still not 
shown. This precursor will not significantly warm the gas, 
however, as the heating occurs where the drag is strong 
and the neutrals decelerate (see Fig. 4. of Smith (1993b)). 

A second problem is inherent to the ion treatment 
in this code: steep ion velocity gradients generate nu- 
merical instabilities (Mac Low & Smith 1997b). Unfor- 
tunately, such steep gradients are exactly what occurs in 
flows set up along the field lines. Hence, at early times, 
up to ~ 6 10 8 s, numerical instability is evident within the 
ion flow. The instability is damped as the gradients di- 
minish. This instability does not occur when, instead of 
discontinuous initial conditions, we use a hyperbolic tan- 
gent function with a width of 1 10 14 cm, and we get nearly 
identical results. 

The tranverse ion motions are found to have a max- 
imum, as was shown to be true for shocks with M a 
> 3-^/2/2, in the steady cold solutions of Smith (1993a). 

We have also set up the quasi-hydrodynamic, or inter- 
mediate Type II boundary and initial jump conditions. In 
intermediate shocks, the post-shock sound speed is sub- 
sonic. As predicted by analytic theory, this flow remained 
J-type and steady. It is clear that a C-type flow should not 
occur since the ions, without the inertia of the magnetic 
field, are now strongly tied to the neutrals. Such Type 
II J-shocks in multiple dimensions could degenerate into 
a fast (switch type) and slow shock combination; exactly 
what is possible depends on the shock configuration and 
applied conditions (see Smith 1993c). We find, however, 
that slow shocks (e.g. Smith 1993c) are also restricted to 
the J-shock variety. In this case, the sub-Alfvenic motions 
tie together the ions and neutrals. 

Type I intermediate shocks are similar to switch shocks 
except the small transverse field reverses direction within 
the shock layer. We find the evolutionary behaviour fol- 
lows that of the switch shock, but at a slower pace (ap- 
proximately 3 times slower for the standard conditions 
with the field initially at 1° to the shock normal.) 

We conclude that a continuous range, from C-type 
to J-type occurs, when the boundary conditions are al- 
tered from switch, via Type I and Type II intermediate, 
to hydrodynamic type. This is emphasized through Fig. 9, 
which is an example of an intermediate shock close to the 
I-II transition border. Here the jump shock remains while 
a partial ion/magnetic precursor develops. Note that here 
also the transverse ion Alfen waves propagate far upstream 
(the complete computational grid is not shown). 

3.5. Neutralisation in transverse shocks 

A class of physical conditions produce another type of 
shock structure in which the ions are not conserved. A 
particularly simple case, in which the ion fraction rc- 
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Fig. 9. The evolution of an intermediate shock structure with 
the field at an angle of 6° to the shock velocity, a compression 
of 57.14, Mach number M = 10, and neutral Alfven velocity 
M a = 5. The parallel field is constant (0.341 mG) while the 
transverse field flips from 0.0362 mG through zero to -1.57 mG. 
The four lines are as described in Fig. 7, and the grid size is 



mains fixed throughout, is considered here. This can oc- 
cur when the recombination time is short, such as for 
molecular ions, with rate coefficients R n <~ 1CP 7 cm 3 s _1 
rather than atomic ions; and when the shocks are wide 
enough for neutralisation reactions to occur (Flower et 
al 1996). Since the shock width is inversely proportional 
to the ion density, neutralisation is determined solely by 
the recombination rate coefficient and Alfven number. 
The ratio of shock dynamical to recombination times is 
~ 100(i?„/l(T 7 cm 3 s-^/Ma. 

We find that neutralised shocks evolve from the im- 
posed J-type initial state directly to C-type. The neutral 
jump section is preceded by a ramp of growing amplitude, 
as shown in Figure 10. The final C-shock flow pattern is 
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now the ramp structure in the neutrals, preceded by the 
rapid ion braking. In this zone of rapid ion braking and 
compression the recombinations will be hard pressed to 
suppress a rise in ion density. This may lead to a narrow 
peak in the ion density in reality. 
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Hence, low-ionisation shocks may never actually reach a 
steady state. 

4. Observables 

Shock waves in a magnetized molecular gas will be con- 
tinuously changing, both through the C-shock instability 
(Wardle 1990, Neufeld & Stone 1997, Mac Low & Smith 
1997a, b) and through interaction with spatial perturba- 
tions in density, ion fraction and magnetic field. So how 
can we recognise the evolving shocks described in this pa- 
per? To calculate the emission properties as a function 
of time, we must first estimate the temperature distribu- 
tion through the shock. In the continuous sections this can 
be done using the 'cool C-shock' approximation (Smith & 
Brand 1990a), provided that the temperature remains low 
so that the thermodynamics and magnetohydrodynamics 
are decoupled. Quantitatively, we require that the temper- 
ature T <C m(H 2 )Ws/k. The low isothermal temperature 
adopted ensures this condition in our numerical computa- 
tions, and strong H 2 cooling often ensures it in real molec- 
ular clouds. 

The jump sections are dealt with separately: the jump 
parameters are transferred into the J-shock code of Smith 
(1994a), the column densities of ions in each excited level 
are calculated, and then added to the column densities de- 
rived from the continuous section in the manner described 
below. 

The temperature profile in the continuous section is 
determined by the local balance of ion-neutral frictional 
heating with molecular cooling (see Smith 1993b). Simpli- 
fied cooling functions are adopted consistent with the lim- 
itations of the calculations already introduced (in particu- 
lar the simplified drag formula). Here we restrict the illus- 
trative results to that of H2 cooling, with the H2 molecule 
in local thermodynamic equilibrium. The cooling function 
is extracted from Smith (1993b): 



Fig. 10. The flow patterns for a neutralised shock transition. 
All parameters are as chosen for the standard transverse shock 
except the ion density is held at its initial value. Note the 
broader shock. The grid size is R g = 200 zones, so the resolution 
is only R = 68 zones in the shock width. The neutral (solid line) 
and ion (dash-dotted line) velocities are shown, as well as the 
momentum error e m (dotted line). 



A(T) = (4.2 lO^ 1 erg s" 1 cm- d )n(H 2 )T 



-.3.3 



(11) 



It is straightforward to show that the temperature is 
then given by T = 2.89niV^ n where Vj„ is the ion-neutral 
streaming speed. 

Column densities of molecules in excited upper levels 
can be computed from 



N,=g 3 N(H 2 )Z(T)eM-T,/T), 



(12) 



The time to reach steady state in a neutralised shock 
is much longer than for the other shock types. Almost 
10 4 years is required for the conditions chosen, compara- 
ble with the duration of the Class stage of a protostar. 
Shock widths are correspondingly larger. Furthermore, the 
age and shock width arc inversely proportional to the ion 
density, as we have checked through further simulations. 



assuming that the rotational levels are in LTE at the tem- 
perature T, where the partition function is 



Z(T) = 0.024 T [1 - cxp(— 6000/T)] \ 



(13) 



and gj are the statistical weights. Although wc assumed 
the rotational levels arc in LTE, we allowed the vibrational 
levels to fall out of LTE using the method described by 
Suttner et al. (1997). Line strengths Ij are then calculated 
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from the column of gas Nj(x) in the transition's upper 
energy level: 

Ij = J (hc/\ J )A 1 N J (x)dx, (14) 

where Aj is the radiative coefficient and Xj the wavelength 
of the transition. 



CO 

C7> 
l_ 
CD 



o 




4x10 



8x10 s 



time (seconds) 



or 
Q 
(J 

cn 
o 



or 
O 



1.2 
1.0 

0.8 
0.6 
0.4 
0.2 
0.0 
-0.2 
1.2 
1.0 

0.8 
0.6 
0.4 
0.2 
0.0 
-0.2 



v s = 30 km s 
HH91A doto 




v s = 40 km s" 1 
HH91A doto 




1x10 4 
Upper Energy, Tj (K) 



2x10 4 



Fig. 11. Line intensities versus time for a 30 kms -1 transverse 
shock model. The shock is assumed to be observed face on (thus 
the columns are calculated along a line of sight parallel to the 
shock velocity). The physical data are given in Table 2. 



Fig. 12. CDR diagrams (see text) for the standard transverse 
field model with the indicated shock speed and evolution times 
of 210 8 s (long dash), 10 9 s (dot-dash), 2 10 9 s (dash), 10 1() s 
(full lines). The data points displayed here were calculated by 
Smith (1994b) from the intensities observed by Gredel et al 
(1992). Only positive line detections are shown here. 



Table 2. Physical data for selected molecular hydrogen emis- 
sion lines, observable in either the K-band (2pim-2.4/im) or by 
the SWS of the Infrared Space Observatory (0-0 transitions). 
Their intensities are plotted in Fig. 11. 
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"Wavelength in /im 
^Excitation temperature in K 

c Einstein A-values for radiative deexcitation in 10~ 7 s 
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The time dependence of eight lines observable cither 
in the K-band or with the SWS of the Infrared Space 
Observatory are shown in Fig. 11, and the physical pa- 
rameters of these lines are given in Table 2. The columns 
of gas are displayed as column density ratios (CDRs) in 
Fig. 12. A CDR for H2 is the column of gas in the en- 
ergy level Tj divided by the column of gas in the v = 1, 
J = 3 level that generates the 1-0 S(l) line, further nor- 
malised by the factor exp(T, /2000K) to remove the strong 
temperature dependence. We thus are comparing column 
densities to those of a slab of molecular gas at a temper- 
ature of 2000 K. These CDR diagrams are an accurate 
means of displaying the H 2 excitation over a broad range 
of energy levels (see also Mac Low & Smith 1997b). A 30 
kms -1 model was chosen for display purposes since this 
produces gas up to a maximum temperature of <~ 2600 K 
within the final steady C-shock. The CDRs for a hotter 
40 kms -1 model are also shown. The following stages are 
recognisable. 

— The structure begins as a fast dissociative jump 
shock with low luminosity in these infrared lines 
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(Fig. 11). The ion-magnetic precursor develops rapidly 
in the early stages, producing a sharp increase in line 
strengths. 

— The jump speed drops below ~ 24 kms -1 at <~ 4 10 8 s. 
Now the molecules survive the jump. The molecules are 
strongly heated in the jump shock producing a high- 
excitation spectrum, as best illustrated by the high 
CDRs during this stage in the high energy levels. Note 
that HH91A possesses such a high excitation spectra, 
as shown in Figure 12. 

— The jump weakens and is unable to excite the 
molecules after <~ 3 10 9 s. The shock excitation and 
line intensities are now controlled by the ion-neutral 
drag and rapidly approach their final values. The final 
excitation depends strongly on the shock parameters. 

We find no significant variations from the above be- 
haviour for the other types of shocks. Switches, oblique 
shocks and intermediate shocks are of course hotter than 
the equivalent transverse shocks. Neutralised shocks pos- 
sess similar excitation signatures as transverse shocks. 

5. Conclusions 

We have simulated the evolution of jump shocks into con- 
tinuous shocks. We set up boundary and initial condi- 
tions appropriate to the C-shock, but separated the two 
flow regimes by a discontinuous diaphragm as in classical 
shock tube experiments. The questions we hope to an- 
swer, however, are more general: how do shocks behave in 
a non-uniform medium? Do steady-state C-shocks form? 
We have indeed found that, in all cases, analytic steady- 
state solutions are approached, given sufficient time. That 
time can in the most extreme cases, however, be compa- 
rable to the lifetimes of embedded protostars. 

Four stages were identified for most shock types. First, 
an ion expansion wave rapidly advances into the upstream 
region. Next, the ion-magnetic precursor forms and moves 
gradually upstream, with a high streaming speed capable 
of heating the molecules. The neutrals develop a weak pre- 
cursor. Then, the neutral jump weakens and disappears. 
Finally, the C-shock runs through minor adjustments and 
reaches the steady state. 

In intermediate Type II shocks (which are quasi- 
hydrodynamic), the full jump shock remains with an ion- 
precursor in advance. Between Type I and Type II, par- 
tial jumps remain. Switch-on shocks develop extremely 
long ion-magnetic precursors due to forward-moving un- 
damped ionic Alfven waves. Switch-off shocks, like all slow 
shocks, are J-shocks. 

We have modelled the time-development of molecu- 
lar hydrogen emission lines in the infrared. If the initial 
jump shock is fast, the molecules do not survive, and lit- 
tle emission occurs. As the precursor develops, the line 
intensities rise rapidly. As the jump shock decreases in 
strength, the molecules begin to survive and emit with 
the high excitation characteristic of a jump shock (e.g. 



with a 1-0 S(l)/2-l S(l) intensity ratio of 3-5). Later, the 
jump is too weak to excite H2 and the C-shock dominates 
the spectrum, producing line ratios that are very sensitive 
to the particular parameters of the shock. 

The evolution time scale is inversely dependent on 
the ion density. For typical molecular cloud parameters, 
it takes between (10 9 /nj) s and (210 10 /ni)s to set up 
the steady-state flow, several times the flow timescale of 
(7 10 8 /ni) s. Neutralised shocks, modelled here by fixing 
the ion number density, are approximately M times wider 
than the transverse equivalent and take the longest to 
evolve. It follows that detectable changes to the emission 
lines occur over periods of at least 60 years for the ion den- 
sity taken here (n^ = 0.1cm -3 ). However, faster evolving 
shocks, due to high ion densities (perhaps in the dens- 
est cloud regions), will show detectable changes faster. 
The time scale for emission-line variations is ~ (2 10 8 /nj) s 
where n^ is the pre-shock ion number density. 

This study provides the background to investiga- 
tions of shock interactions, C-type jet flows, and multi- 
dimensional studies such as simulations of the Wardlc in- 
stability (Mac Low & Smith 1997a, b). We hope to add ad- 
ditional physics to this version of ZEUS, including chem- 
ical reactions, molecular dissociation and streaming ioni- 
sation. 

Acknowledgments: MDS thanks the DFG for financial 
support. 
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